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Abstract. We derive the first-passage-time statistics of a Brownian motion driven 
by an exponential time-dependent drift up to a threshold. This process corresponds to 
the signal integration in a simple neuronal model supplemented with an adaptation- 
like current and reaching the threshold for the first time represents the condition for 
declaring a spike. Based on the backward Fokker-Planck formulation, we consider the 
survival probability of this process in a domain restricted by an absorbent boundary. 
The solution is given as an expansion in terms of the intensity of the time-dependent 
drift, which results in an infinite set of recurrence equations. We explicitly obtain 
the complete solution by solving each term in the expansion in a recursive scheme. 
From the survival probability, we evaluate the first-passage-time statistics, which itself 
preserves the series structure. We then compare theoretical results with data extracted 
from numerical simulations of the associated dynamical system, and show that the 
analytical description is appropriate whenever the series is truncated in an adequate 
order. 



FPT statistics of a Brownian motion driven by an exponential temporal drift 2 
1. Introduction 

The statistical analysis of a system is essential when fluctuations contribute to its 
dynamics [DEJE]. ^ n different areas, beyond the importance of the statistical description 
of the system state and its evolution, the main variable of interest is the time at 
which this state reaches a certain region for the first time [I], constituting the so- 
called first-passage-time (FPT) problem. For example, in a diffusion-controlled reaction 
a particle performs a random walk until it makes contact with a reactant or a trap 
giving rise to the reaction [3]. Generally, in a FPT problem, the system is able to 
evolve according to a given dynamics in a confined region, limited by one or more 
absorbing boundaries. Even for simple autonomous systems, the FPT problem can be 
analytically difficult to solve. For example, for the Ornstein-Uhlenbeck process with 
a fixed positive absorbing boundary, the FPT solution (representing the FPT density 
function) is relatively easy to be found in the Laplace domain, but its inverse transform 
is not explicitly available [5]. A notable exception is the FPT problem for a Brownian 
motion (Wiener process), where different methods are easy to be applied to solve it 
[U [21 IH El E]- A greater complexity is found when a time-inhomogeneous process 
defines the system dynamics. In this case, analytical methods are formally given within 
different approaches [21 El El E], but exact as well as approximate explicit results 
are scarce in the literature and probably difficult to obtain. Among the different ways 
to introduce a time inhomogeneity into the system (for example, temporally varying 
absorbing boundaries [HI EE], time-dependent drift and diffusion coefficients [UJ, etc), we 
focus on a particular drift coefficient evolving externally in time. The process analysed 
here is driven by an exponential time-dependent drift (this case, in turn, can be mapped 
into a variable threshold [9]), which naturally arises in neuroscience when modelling 
adapting neurons [121 H31 [HI [151 EE El • m this context, the membrane potential (state 
variable) of a perfect integrate-and-fire neuron model is driven during its subthreshold 
evolution by an external current, composed of a constant deterministic current plus fast 
fluctuations [21 El El [321 EES] , as well as an intrinsic temporally decaying current [T71 [20] ■ 
This system corresponds exactly to the case analysed here, and the FPT represents the 
production of a spike. 

The interest in the FPT problem for a Brownian particle in a time-inhomogeneous 
setup started in the 1990s, when different systems driven by periodically modulated 
drifts were studied within the context of the stochastic resonance phenomenon [21]. In 
particular, Bulsara et al applied the method of images to a Wiener process driven by a 
sinusoidal temporal drift in the presence of an absorbing boundary [22J , a procedure of 
limited validity [231 EH E] • Later, other threshold processes under analogous conditions 
were theoretically analysed with different approximation methods [211 [23 [2S1 [2Z] • Even 
when appealing, single sinusoidal temporal drifts do not represent a general case. For 
arbitrary time-dependent drifts, the simplest procedures are a quasiadiabatic reduction 
[28J, a quasistatic description [29J, or small amplitude approximations. However, for 
rapidly varying arbitrary fields, the time-dependent structure of the problem cannot be 
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simplified. Up to our knowledge, the first attempt to include an arbitrary temporal 
drift (without spatial dependence) in a general framework was made in [30J, where the 
author proposed a method to describe the first-order correction to the moments of the 
FPT density, in a perturbation scheme. 

The preceding studies describe approximately the FPT problem of a given 
continuous stochastic process in a restricted domain for particular or general temporal 
drifts. In general, explicit exact results for time-inhomogeneous systems are infrequent. 
Notably and as an exception, in [IJJj the authors derive the FPT density function of 
a Wiener process, in the presence of an absorbing boundary, where both the drift 
and the diffusion coefficients are varying temporally and in proportion to each other. 
In particular, when proportionality is satisfied, the Fokker-Planck equation ruling the 
evolution of the transition probability between the states at two times can be time- 
rescaled in order to resemble the simpler constant coefficients case, where the exact 
solution is known. However, the restriction in the coefficients proportionality limits its 
applicability to our case. For the FPT problem we are interested in, we have previously 
proposed a series solution in terms of the intensity of the time-dependent drift [31] (see 
also [221 EH] for analogous series solutions, but focusing on the perturbation regime). 
However, in that work only the first terms in the expansion were explicitly given and 
higher order terms were just outlined. 

In this work, based on the structure of the equations obtained in [31] for the survival 
probability, we explicitly obtain all order functions. In particular, we find the nth- 
term in a recursive scheme and prove by induction the complete mathematical solution. 
From the survival probability, it is straightforward to derive the FPT statistics, which 
maintains the series structure. Since obtaining all order functions (existence) does not 
imply the convergence of the series, in the second part of the work, we analyse how 
does the expansion behave in comparison with results obtained from simulations, as we 
truncate the series at a finite order. 

2. Theoretical framework 

In this section, we set the system under analysis, the formalism we use to study the 
survival probability and the FPT density function, and derive the complete solution. 

2.1. The system 

The dynamics of the system is governed by the Langevin equation 

^ = M+-e-^+e(t), (1) 
dt r d 

where x is the state variable (e.g. particle position, membrane potential, etc), t is the 
time, \i is the constant (positive) component of the drift, e and r d characterize the 
intensity and the time constant of the exponential time-dependent drift, respectively, to 
sets the initial time and £(t) is a Gaussian white noise with (constant) squared intensity 
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D [<£(*)) = and mm) = 2D *(t ~ t')}. 

Following the derivation we have made in [31], the probability that the particle 
remains in the domain x < Xthr at time t', given the (variable) initial condition x at 
(variable) time t, F(t'\x,t), evolves according to the backward FP equation 



dF(t'\x,t) 
dt 



H H e 



-(*-*0)/Td 



dF(t'\x,t) d 2 F(t'\x,t) 



dx 



D 



dx 2 



(2) 



In Eq. ((21), is a parameter accounting for the present time. By making the 
substitution r = t' — t and renaming the probability as F(x,r;t'), Eq. ^ can be 
written as 

d 2 F(x,T;t') 



dF(x,r;t') 
dr 



e -(t'-t )/T d g-r/Td 



Td 



dF(x,r;t') 
dx 



+ D 



dx 2 



-,(3) 



which represents the equation to be solved. The system is completed by specifying the 
initial and boundary conditions [3T], which are 



F(x,r = 0;t') = 
F(x = x thr ,r;t') = 0. 



1 for X < Xthr, 
for X > Xthr, 



(4) 
(5) 



Equation Q indicates that the survival at initial time is certain for a particle 
located in the domain of interest, whereas Eq. ^ establishes that the particle is not 
allowed to be in x > x t hr and, therefore, x = x t hr is an absorbent boundary. 

By proposing a solution as an expansion in e, 



F(x, t- t') = F (x, r; t') + e F x {x, r; t') + e 2 F 2 (x, r; t') + . . . = £ e" F n (x, r; 0, (6) 
Equation ^ reads 



n=0 



dF 8F d 2 F 
(j,— D 



Or 



dx dx 2 
dF n 



n=l 



_ J^LJl _ — e -(t'-t )/r d e r/r d 

dr dx Td dx 



dF n ^ D d 2 F n 



dx" 1 



0. (7) 



Since we expect that all functions F n do not depend on e, Eq. ([6]), the expressions 
between brackets should be identically 0. Obviously, this hypothesis is true if we are 
able to find F n (x,r;t'). Under this condition, the complete solution for F(x,r;t') is 
given by the system of equations 



dF 




dF 


- D 


d 2 F 


dr 


- V 


dx 


dx 2 


dF n 




dF n 


- D 


d 2 F n 


dr 


- V 


dx 


dx 2 



jL e -(f-*o)/T d e r/r d dF < 
Td 



n-l 



dx 



for n > 1. 



(8) 
(9) 
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Consistently with our previous assumption, given the arbitrariness of e, the non- 
homogeneous initial condition, F(x, r — 0; t r ) — 1 for x < x thr , should be exclusively 
imposed to the zeroth-order function Fq(x, t = 0; t'). In detail, initial conditions are 



F (x,r = 0;t') 



1 if x < x thl , 

if X > Xthr, 

F n (x,r = 0;t') =0 for n > 1. 



Completing the description, the boundary condition reads F n (x 
for all n. 



(10) 
(11) 

£thr, r; t') = 0, 



2.2. Survival probability from the backward state 

To obtain exactly the survival probability at time t' from the backward state we have 
to solve all terms involved in the expansion given by Eq. ([6]). In particular, each 
term satisfies a certain equation, Eq. ([8j or ([9]), with appropriate conditions. Due 
to the different mathematical structure, we focus on the zeroth-order term, Fq(x, r;t'), 
separately from all other superior terms, F n (x, r; t') for n > 0. 



2.2.1. Zeroth-order term. This term corresponds to the survival probability at time 
t! of a Brownian particle (initially) located in x at time t, when the system is driven 
exclusively by a constant positive drift fi (in our system, this is obtained with e = 0). 
According to the preceding derivation, F (x,r;t') satisfies Eq. (|sj> with the conditions 
given by Eq. (10) and F (x t hr, t; f) = 0. Since F (x,r;t') = for x > x t hr, we focus 



exclusively om< x t hv'i m this case, by Laplace transforming Eq. rt8j), we obtain 



s F L (x) - n 



&Fk(x) D d 2 F n L ( 



dx 



dx 2 



(12) 



30 e ST F (x,T;t') dr is the Laplace transform of F (x,T;t') in 



where Fq(x, s;t') = / ( 
the variable r. Since t' and s act as parameters in Eq. (12), we have simplified the 
notation to Fq(x). In Laplace domain, the boundary condition simply transforms to 

i^thr) = 0. 



The solution to Eq. (|12|), with the preceding condition and taking into account that 

), is 

(^thr — x ) 



Fq{x) remains bounded asn> — oo, is 



Fn(x) 



1 1 

exp 



2D 



(13) 



The inverse Laplace transform of Fq (x, s; t') can be explicitly computed and reads 



- erfc 
2 



exp 



X 



2VDt 

(X th r - X) (i 



D 



H T 

2VD 
erfc 



(^thr 



X 



2VDt 



+ 



(14) 



where erfc(x) is the complementary error function. 
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2.2.2. Higher order terms. The equation governing the dynamics of the nth-order 
function is given by 



dF n (x,r;t') dF n (x,r;t') d 2 F n (x,r;t') _ 1_ _ (t /_ to)/rd 



dr 



dx 



dx' 



Td 



d_ 
dx 



e r/Td F n ~i(x, t; if)] , (15) 



and the boundary and initial conditions are F n (x t h T , t; t') = and F n (x, r = 0; t') = 0, 
respectively. 

This equation can be solved via a Laplace transformation in the variable r, which 
reads 



sF^(x)-fi 



dF^(x) 
dx 



-D 



d 2 FHx) 1 



-(t'-t )/T d 



dx 2 



dx 



s-l/r d 



(16) 



where F^(x,s;t') is the Laplace transform of F n (x,r;t'), and its notation has been 
simplified to F£(x) as in the previous case. Due to the exponential pre- factor, the 
Laplace transform of the forcing term has to be evaluated in the shifted variable s — l/r d . 
Again, the boundary condition is simply F^(x thr ) = 0. 
Next, we prove that the nth-order function is 



-n(t'-t Q )/T d _ 



[pi-^fi 2 + 4D(s-n/r d )] 
2D(s-n/T d ) 



x J2 a nA s ) ex P 

i=0 



( x thr ~ x ) 

2D 



[^-J^ + AD(s-i/T d )} 



(17) 



where the n + 1 coefficients weighting each of the exponential terms appearing in the 
solution of the nth-order function, a n ^(s) (i — 0, . . . , n), are given by 



f , A a„-i,i-i(5 - l/r d ) \» ~ V^ 2 + AD ( S ~ *Ad)] 
1=1 1 



2D 



a n .i(s) 



^(s-l/rd) [fi- Jfi 2 + 4,D(s-i/T d 



2D 



for i = 1, 



(19) 



which build a recursive solution. In particular, from Eq. (19) it is easy to check that 

To demonstrate this solution, we will set F r f_ 1 (x) according to the preceding 
proposition and prove that the following order function satisfies the same structure, 



Eq. (17). In doing so, we will derive explicitly the recursive scheme given by Eqs. (18) 



and (19). The demonstration will be completed by showing that the first-order term, 
F[(x), belongs to the family of functions defined by Eq. (17) (i.e. we prove the 
proposition by mathematical induction). 

Given that F^_ ± (x) is expressed according to Eq. (17), 
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^-i(z)=e- 



-{n-l){t'-t )/T d 



{H- Jfi 2 + AD[s-(n-l)/r d ]} 



n-l 

x 0>n-l,i( s ) ex P 

8=0 



2D[s-{n-i)r d ] 

{%thr — X 



2D 



[fi-Jfi 2 + AD(s-t/r d )} 



(20) 



the next order function, F^(x), is given as the solution of Eq. (16) with an explicit 
forcing term, 



dFMx) „ d 2 F„ L (aO 



dx 2 



— e 



-n(t'-t )/T d 



[/J - tJ/j, 2 + 4D(s -n/r d ) 
2D r d (s - n/r d ) 



n-l 



X a n-l,i( S ~ l / T d) 



8=0 



x exp 



(Xthr ~ 

2D 



{/i- v //i 2 + 4D[ g -(^ + l)/r d ]} 
2D 

{/i - v /^ + 4jD[s _ (z + i )/rd ]|| . 



The solution to the homogeneous part of this equation reads 



^homfc) = Cl eX P 



-/i + y//i 2 + ADs 

2D 



x \ + C 2 exp 



-/i - V/i 2 + ADs 
2D 



whereas it is easy to check that a particular solution is 



(21) 



x , (22) 



[/i- ^ 2 + AD(s - n/r d ] 
2D(s - n/r d ) 



F L (t\ = _p- n (*'-*0)/T d 

- 1 n^partl" / c 



x g a n _ M (s - l/r d ) {(* - yV + 4D[s - (z + l)/r d ]} 



8=0 



z' + l 



2D 



x exp { (Xth 2D X) ^-^ 2 + 4D[ S -(z + l)/r d ] 



(23) 



The general solution is obtained from the combination of Eqs. (22) and (23), 



F„ (x) = F^ hom (x) + F^ paTt (x), and it is valid for Re(s) > n/r d . Since F^(x) is bounded 
as x — )■ — oo, C 2 is 0; at the same time, the boundary condition, F% (xthr) — 0, builds 



d = e 



-n(t' -to)/Td 



[u - J/i 2 + 4D(s - n/r d )} ffi- V/i 2 + ADs 



2D(s - n/r d ) 



exp 



2D 



n ^ a n ^{s - l/r d ) {/i - V/" 2 + *D[s - (i + l)/r d ]} 



z' + l 



2D 



(24) 



Therefore, the nth-order function is given by 
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-n(t'-t )/T d 



xj<Vo( s J exp 

n-1 

+ a n,i+i(s) exp 

i=0 



\p - y/i 2 + 4D(s - n/r d ; 
2D(s-n/r d ) 

(^thr ~~ x 



2D 

j (%thr ~ x 

I 2D 



{u- yW^s 



(i + l)/r d ] , (25) 



where the coefficients appearing in Eq. (25) are 



, , ^ a n ^(s - l/r d ) {H ~ y/fi 2 + W[s-(i + l)/r d ]} 

a n n (s = > 

'° V ' i + 1 2D 



i + 1 



2D 



(26) 



On- M (a ~ lAd) {/* ~ yV 2 + 4 ^ - (» + l)/r d ]} . 
a n>i+ i(sj = y - — , z = 0,...,n-l. (27) 



By shifting the index z in the sum symbol, it is easy to check that Eq. (25) is equal 



to Eq. (17), and each of the coefficients, Eq. (26) or Eq. (27), is given by Eq. (18) or 



Eq. (19), respectively. 



The proof is completed by showing that the first-order function, given as the 
solution to Eq. (16) for n = 1 and F^Xthr) = 0> * s P ar ^ °f the family of functions 



described by Eq. (17). As shown in [31], this solution reads 



[fi-J^ + AD(s-l/r d ) 



2D(s - l/r d ) I 
+ exp | {Xth ^ D X) [n - y/f + W(3-l/T d )] 



exp |(xt^ [ ^_ v r^^ ] 



(2* 



which can be easily checked satisfying Eq. (17). Furthermore, from this solution we can 



observe that the coefficients a n ^{s) (n — 2, . . . , oo and % — 0, . . . , n) are recursively built 
from aifi(s) = —1 and Oi,i(s) = 1. 

Even when not explicitly available, the nth-order function in the temporal domain, 



F n (x, t; t'), is given by the inverse Laplace transform of Eq. (17), which reads 



F n (x,r;t') = e 



-n(t'-t )/T d 



1 r +jo ° rr&- J/-i 2 + ±D(s-n/T d )) 



2it] Jo--joo 2D(s — n/r d ) 

x lt a nA s ) exp | ^ tbr 2D X \ n-^t> 2 + 4-P(s - i/r d )] j ds, 



(29) 



where j represents the imaginary unit and the region of convergence of the integrand 
requires that o > n/r d . From the substitutions z = s — n/r d for the integration variable 
and k = n — i for the index of the sum, we obtain 
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F n (x,r;t') 



e -n(t'~t )/T d e nr/T d 



1 r^°° ZT [fi - y/n 2 + ADz 

27TJ Ja z -]oo 2Dz 



\ (x t hr ~ 



x) 



(30) 



where now, o~ 2 > 0. In Eq. (30), we have defined new coefficients for the exponential 
terms appearing in the sum symbol, b n> k(z) = a n ^ n _ k (z + n/r d ). With this definition, 
the recursive structure is 



n-k 2D 

n-l 

b n ,n(z) = ~ b n,k( Z )' 

starting from 61,0(2) = 1 and 61,1(2) = — 1. 



, for k = 0, . 



.,n-l, 



(31) 
(32) 



2.3. Survival probability 

The survival probability of the particle at time if arises when we impose the initial 
state to the backward state, x = xo at time t = t . As in the previous subsection, we 
discriminate between the zeroth-order term from all superior order functions. 



2.3.1. Zeroth-order term. This term is given by imposing the initial state in Eq. (14) 
and reads 



FJt) 



- 1 , rOrthr - x ) /i f?" 
1 erf c , 4 / — 



exp 



(Xthr - go) jl 

D 



erfc 



(X t hr ~ Xq) 

2\fm 



fi [T 

2VD 



(33) 



where now, r = t' — 1 is the actual time difference (time elapsed from the initial time 
to to the present time if). Note that we have eliminated the dependence on t' in the 
notation for F (t), since it only appears in the combination given by r. 
The Laplace transform of Fq(t) in the variable r is given by 



1 1 J(Xthr-X ) r / 2 j , n 1 



(34) 



which is one of the quantities of interest for the assessment of the FPT density function. 
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2.3.2. Higher order terms. As in the zeroth-order, these terms are obtained from the 
evaluation of the initial state in the corresponding expression for the survival probability 



from the backward state, Eq. (30), which yields 



27TJ Ja z -)oo 2Dz 

\ 2D 



x J2 Kk(z) exp { v ™ n u > - V/i 2 + AD(z + k/r d )} \ dz, (35) 

fc=0 



where r = if — to is the actual time difference and the dependence on if appears only 
through r. 

The Laplace transform of F n {r) in the variable r is readily obtained, and reads 

~ L _ [fi-V^ + 4Ds} 
n K ) ~ 2Ds 

x jrb rhk (s) exp | {Xthl 2 ~ Xo) [ft - ^ + AD(s + k/r d )] | , (36) 



where the coefficients b nj k(s) are given by Eqs. (31) and (32) in the variable s. 



2.4- First-passage-time statistics 

The probability that a Brownian particle driven by an exponential time-dependent drift 
(superimposed to a linear field) remains in the domain x < x t hr a t time t', having 
started at time to in the position xq, is given by the survival probability calculated in 
the previous subsection, F(r). It was demonstrated that this probability can be written 
as a series 

oo 

F(r) = Y,e n F n (r), (37) 
n=0 

where each term depends exclusively on the variable r = if — to. The different order 



terms are obtained from the functions explicitly found in Eqs. (34) and (36), in the 
Laplace domain. 

For a positive fi, the particle will cross the level x = Xthr for the first time at time 
T (and will be absorbed), and this random variable represents the FPT. Given the 
survival probability at time r (time elapsed from time to) ; the FPT for this particle 
satisfies T > r (it is absorbed at a posterior time); therefore, the survival probability 
represents F(r) = Prob(T > r) and the cumulative distribution function for the FPT, 
<3>(r), is given by $(r) = 1 — F(t). Consequently, the density function for the FPT, 
0(r), is 
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Since F(r) is given as a series solution, Eq. (37), the density function 0(r) can also 
be expressed as a series, 

oo 

= 5> B ^(r), (39) 

n=0 

where the functions <p n (r) are 

dF n (r) 



dr ■ 

In Laplace domain, the preceding series is expressed as 



(40) 



oo 

IL I „\ \ n ~lL 



> S 



£^#00, (41) 

n=0 



where [31] 



^( S ) = 1- S F L ( S ), (42) 
5£(s) = -s P£(s), forn>l. (43) 



Replacing the results we have obtained in the previous subsection, Eqs. (34) and 



(36), these functions explicitly read 



exp I (Xth ^ Xo) [/i - (44) 
[/i - V/i 2 + ADs 



x 



2D 

f: 6 Bifc (s) exp | {Xthr 2D Xo) \ji - y/n* + W(s + k/T A )]} , forn > 1, (45) 



where the coefficients b nt k(s) are given by Eqs. (31) and (32) in the variable s, and the 
recursive structure starts from &i,o(s) = 1 and &i,i(s) = — 1. In order to exemplify this 
recursive construction, in table [T] we show the coefficients associated with the exponential 
terms, b n ^(s), up to the fourth order. 

3. Comparison to numerical simulations 

To illustrate the solution we have obtained for the FPT problem of a Brownian particle 
driven by an exponential time-dependent drift, in this section we compare the analytical 
results with data extracted from numerical simulations. Samples of the FPT distribution 
are collected from the times at which a Brownian particle, evolving according to the 
Langevin equation given by Eq. and starting from xq, arrives at the threshold x t hr 
for the first time. As we have shown in [31], for small intensities of the time-dependent 
drift, e, the system corresponds to a perturbation scenario, and the first-order solution, 
4>(t) = 4>o(t) + e (f>i (r), properly describes the FPT statistics. In this section, we extend 
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that comparison beyond the linear regime, for large values of e. Obviously, in this case 
we need to include higher order terms in the series given by Eq. (37). This comparison 
is not merely illustrative; since the characterization of the series convergence remains 
elusive to us, we resort to a numerical test case to demonstrate the usefulness of the 
series solution. 

Considering e = 0, the dynamics defined by Eq. ([I]) has no intrinsic timescale 
and, therefore, time units can be normalized by (scthr — x o)/t L (i- e - the external 
parameter /i defines the escape rate). Additionally, we consider the non-dimensional 
form of this equation, obtained by setting x/(x thr ~ x o) ~ > x i e /( x thr — x o) e > an d 
D/[(xthx — x o)fj] — > D. The preceding procedures are equivalent to set /i = 1 and 
x thr— x o = 1 in the system described by Eq. ([I]). Given our interest in neural adaptation, 
the remaining parameters will be defined from typical values in adapting neurons. For 
the system without the time-dependent drift, e = 0, the mean FPT is (r) = 1; in the 
context we focus on, a proper scale for r d is about 10 times this value [16], = 10. 
Since the squared noise intensity strongly influences the dispersion of the interspike 
interval distribution (FPT statistics), its value is selected to produce typical histograms 
obtained in experiments pj, D = 0.01. The intensity of the time-dependent drift weights 
the influence of an adaptation current in the intrinsic FPT distribution (in particular, 
in the firing rate) and constitutes a negative feedback to the subthreshold integration 
[32] (e < 0); given that our interest here is to analyse the behavior of the series solution, 
this parameter will be used to set different regimes beyond the linear case. 

In Figs, ([lja) and (|T|b) we show the FPT density function constructed from 
numerical data, for different intensities of the time-dependent exponential drift, e. As 
expected, for positive (negative) intensities, as the strength of the time-dependent 
drift increases in magnitude, the threshold is reached at earlier (later) times and, 
consequently, the FPT distribution shifts towards lower (larger) values. In Figs. (|T|c) 
and ([l}d), the different distributions are separately compared with analytical results. 
These results are based on the truncated series solution [see Eq. (37)], 0(r) = 
J2n=o e " ^n( T ), where N is selected to reproduce numerical data. As shown in Figs. ([l]c) 
and (jTJd) (top panel), the FPT distribution for e = ±0.1 is precisely described by the 
linear expansion (perturbation regime), <p(r) = 4>q{t) ± e <pi(r). However, the proper 
description of the numerical distributions for higher intensities requires the addition of 
higher order terms. For example, as shown in Figs. (Jljc) and (jTjd) , FPT distributions 
for e = ±0.5, e = ±1.0, and e = ±2.0 are described with N = 2, N = 4, and N = 9, 
respectively (middle-top, middle-bottom, and bottom panels, respectively). 

Except for the zeroth-order term, which can be explicitly computed in the temporal 
domain via the inverse Laplace transformation of Eq. (44) and known as the inverse 
Gaussian distribution [BJ, 

MT> = CXP I Wr ) ■ (46) 
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all superior order functions, 4> n (^) f° r n > 1 ; require the numerical (inverse Laplace) 
transform of Eq. (45). In Fig. ^ we show these functions up to the eighth order, for the 
parameters defined in Fig. 0. It is worthwhile to note that these functions decrease in 
amplitude as the order increases (see ^/-scales), which is indicative of the convergence of 
the series (but not conclusive). An additional point to take into account in this analysis 
is that the numerical transform introduces an error which limits the reliability of the 
results. In particular, for the test case used here, the numerical inversion of the functions 
beyond the tenth order is inaccurate and, therefore, the comparison between analytical 
and numerical results is restricted to e ~ ±2.0 [Figs, (jljc) and (JTjd)] . 

This numerical inaccuracy can be circumvented if we analyse properties that can 
be obtained directly from the Laplace transform of the FPT density function, <f) L (s); for 
example, its moments read 



<T') = ^(T)T*dr = (-l)*^g£>j^. (47) 



It is easy to check that, due to the linear nature, Eq. (47) adopts a series structure 
when L (s) is replaced by Eq. (41). Explicitly, by defining 

= tL„- («) 



Equation (47) results in 

oo 

(r fc ) = £e"(rV (49) 

n=0 

In Fig. ^ we show a comparison between numerical and analytical results for the 
first four moments as a function of the intensity of the time-dependent drift. Parameters 
are defined as those corresponding to the test case analysed previously. This means that 
all (T k ),j >n are certain scalars and Eq. (49) represents a polynomial in e. It is interesting 
to note that the order of the truncated polynomial, N, necessary to describe the nu- 
merical results increases as the exponent of the moment does. Alternatively, for a given 
order, the analytical expressions for the lowest moments remain valid in a larger range 
of |e|. Also, it is interesting to point out that Eq. (49) implies a different behavior for 
positive or negative values of e. As shown in Fig. ([3]), for e < the convergence of the 
analytical expression is smooth as the order of the polynomial increases, whereas for 
e > the convergence exhibits an alternating character. 

From the moments of the FPT density function we can obtain other proper- 
ties; in particular, its cumulants. In this case, the expressions relating both proper- 
ties should be developed, and a series is obtained by grouping together equal order 
terms. For example, the second cumulant is, up to the first order in e, ((r — (t)) 2 ) = 
(r 2 )^ - (r) 2 J + e [(r 2 )^ - 2(r), o (r) 01 ] + (9(e 2 ). 
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4. Concluding remarks and discussion 

We have studied the FPT statistics resulting from the biased diffusion of a Brownian 
particle up to a threshold x t hr, when the constant drift is supplemented with an 
exponential time-dependent component [see Eq. ([!])]. In a previous work [31], we 
analysed this time-inhomogeneous system in the backward FP formalism, derived the 
diffusion equation governing the evolution of the survival probability from the backward 
state, Eq. ([3]), and proposed a solution as a series in terms of the intensity of the time- 
dependent drift, Eq. (J6|. In that work we focused on a perturbation regime and explicitly 
solved the expansion up to the first-order terms. In this work, we have extended these 
results by explicitly computing all superior order functions in a recursive scheme [see Eq. 



(30)]. The survival probability (with the initial state imposed) and the FPT statistics 



are easily derived from this solution and preserve the series structure; for completeness, 



their superior order terms are also explicitly given [see Eqs. (36) and (45), for the 
corresponding expressions in the Laplace domain]. In the second part of this work we 
have defined a test case in order to assess the usefulness of the series solution. Analytical 
and numerical results are compared for different intensities of the time-dependent drift 
(beyond the perturbation regime), and a remarkable agreement is found for each case 
whenever the series is truncated in an adequate order. 

The problem we have analysed provides the intrinsic statistics of the events defined 
by an adapting neuron (interspike intervals). In this case, the system state corresponds 
to the membrane potential and the exponential time-dependent drift resembles a 
specific ionic current that decays during the subthreshold integration. This kind of 
currents supports a widely observed phenomenon in neurons, known as spike-frequency 
adaptation (SFA), when the initial state of the current (in the present framework, 
proportional to e) is properly coupled with the spiking history [32] . Particularly, they 
are restricted to be negative (e < 0), providing a feedback to the neuron that lengths the 
interspike interval (FPT). In this work, we have focused on the statistics describing a 
single interspike interval for a given initial current [i.e. the FPT statistics analysed here 
corresponds to a conditional distribution, 0(r|e), in the history-dependent spike train]. 
As shown in [22], the analysis of the successive events in a neuron exhibiting SFA can 
be performed with a hidden Markov model. In this case, the conditional distribution is 
essential to study the spike train properties and its explicit assessment has motivated 
the contribution made in this study. 
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Table 1. Coefficients b n ,k(s) weighting the exponential terms that compose the 
solution to the survival probability and the FPT density function up to the fourth 
order. For the sake of clarity, we have defined the auxiliary coefficients c n (s) = 
[/i - ^^+4D( S + n/r d )]/(2D). 



Order 


Coefficients 




n = 1 


&i,o(*) = 1 






6i,i(s) = -l 




n = 2 


M*) = -[&i,o(«)/2]co(s) = 


-|c (s) 




62,1 (*) = -[6i,i(«)/l]ci(«) - 


Ci(«) 




6 2 , 2 («) = -&2,o(s) - 62,1 (s) = 


= i Co (s)-ci(s) 


n = 3 


6 3 ,o(s) = -[6 2 ,o(s)/3]c (s) = 






63,1 (*) = («)/2]ci(») = 


-|[ci(«)] 2 




6 3;2 (s) = -[6 2 , 2 (s)/l]c 2 (s) = 


[-|co(s) + ci(s)]c 2 (s) 




63,3(5) = -63,0 (s) - 63,1(5) - 


- 6 3 , 2 (s) = -i[c (s)] 2 + i[ci( S )] 2 + [§<*(*) - Cl ( S )]c 2 ( S ) 



n = 4 6 4! o(s) = -[6 3 , (s)/4]c (s) = -^[c (s 
6 4 ,i(s) = -[6 3 ,i(s)/3]c 1 ( S ) = i[c 1 ( S )] 



3 



2 



64,2(5) = -[6 3 , 2 (s)/2]c 2 ( S ) - \[\c {s) - Cl (s)][c 2 ( S )] 
64,3(5) - -[6 3i3 ( s )/l]c 3 ( s ) = {i[c ( S )] 2 - i[ Cl ( S )] 2 + [-|c ( S ) + Cl ( S )]c 2 ( S )}c 3 ( S ) 
64,4(5) = -64.0(5) - 6 4i i(s) - 6 4 , 2 (s) - 64,3(3) 

= ^M*)] 3 - \[ Cl { S )f + §hf <*(*) + Cl (5)][c 2 ( S )] 2 
+{-\[c^)] 2 + |[ C i(5)] 2 + [|c (s) - C1 (5)]C 2 ( S )}C 3 (5) 
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First-passage-time First-passage-time First-passage-time 



Figure 1. Comparison between the series solution for the FPT density function 
and numerical results, (a) and (b) First-passage-time distributions obtained from the 
numerical simulation of Eq. ([IJ, for different (a) positive and (b) negative intensities 
of the time-dependent drift e (colored histograms), (c) Each histogram shown in (a) 
is properly described by the series solution, Eq. (37), shown as a thin yellow line, (d) 
Equivalent comparison between theoretical results and the histograms shown in (b), 
for negative intensities. In each case, the series is truncated in an adequate order, N: 
4>{t) — J2n=o en ^n( r )- As the value of e increases in magnitude, the order used to 
represent the theoretical result increases as well. In particular, N — 1 for e = ±0.1 
(top), N = 2 for e = ±0.5 (middle-top), N = 4 for e = ±1.0 (middle-bottom), and 
N = 9 for e = ±2.0 (bottom). Parameters: fj, = 1, x t hr — xq = 1, D = 0.01, and 

Td = 10. 
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Figure 2. Functions <t> n (T) used to construct the theoretical description of the FPT 
statistics for the test case defined in Fig. ([I]). Note that the y-scale varies from panel 
to panel and, particularly, decreases as the order becomes higher. 
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Figure 3. Comparison between numerical and theoretical results for the moments of 
the FPT distribution, as a function of the intensity of the time-dependent drift e. In 
all cases, symbols represent averages of numerical results, whereas lines correspond to 
Eq. (49) truncated at the order indicated. Parameters are defined as in Fig. (JlJ. 



